using DynamicPolynomials
using SumOfSquares
Sum-of-Squares programs are usually solved by SemiDefinite Programming solvers (SDPs). These programs can be represented into two different formats: Either the standard conic form, also known as kernel form: $$ \begin{aligned} \min\limits_{Q \in \mathbb{S}_n} & \langle C, Q \rangle\\ \text{subject to:} & \langle A_i, Q \rangle = b_i, \quad i=1,2,\ldots,m\\ & Q \succeq 0, \end{aligned} $$ or the geometric conic form, also known as image form: $$ \begin{aligned} \max\limits_{y \in \mathbb{R}^m} & \langle b, y \rangle\\ \text{subject to:} & C \succeq \sum_{i=1}^m A_i y_i\\ & y\ \mathsf{free}, \end{aligned} $$
In this tutorial, we investigate in which of these two forms a Sum-of-Squares constraint should be written into. Consider the simple example of trying to determine whether the following univariate polynomial is a Sum-of-Squares:
import SCS
@polyvar x
p = (x + 1)^2 * (x + 2)^2
model_scs = Model(SCS.Optimizer)
con_ref = @constraint(model_scs, p in SOSCone())
optimize!(model_scs)
------------------------------------------------------------------ SCS v3.2.4 - Splitting Conic Solver (c) Brendan O'Donoghue, Stanford University, 2012 ------------------------------------------------------------------ problem: variables n: 6, constraints m: 11 cones: z: primal zero / dual free vars: 5 s: psd vars: 6, ssize: 1 settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07 alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1 max_iters: 100000, normalize: 1, rho_x: 1.00e-06 acceleration_lookback: 10, acceleration_interval: 10 compiled with openmp parallelization enabled lin-sys: sparse-direct-amd-qdldl nnz(A): 12, nnz(P): 0 ------------------------------------------------------------------ iter | pri res | dua res | gap | obj | scale | time (s) ------------------------------------------------------------------ 0| 1.30e+01 1.87e+00 2.95e+01 1.47e+01 1.00e-01 2.84e-04 50| 6.01e-07 2.84e-07 1.95e-06 -9.76e-07 1.00e-01 4.45e-04 ------------------------------------------------------------------ status: solved timings: total: 4.46e-04s = setup: 1.14e-04s + solve: 3.32e-04s lin-sys: 1.92e-05s, cones: 1.84e-04s, accel: 4.05e-06s ------------------------------------------------------------------ objective = -0.000001 ------------------------------------------------------------------
As we can see in the log, SCS reports 6
variables and 11
constraints.
We can also choose to dualize the problem before it is
passed to SCS as follows:
using Dualization
model_dual_scs = Model(dual_optimizer(SCS.Optimizer))
@objective(model_dual_scs, Max, 0.0)
con_ref = @constraint(model_dual_scs, p in SOSCone())
optimize!(model_dual_scs)
------------------------------------------------------------------ SCS v3.2.4 - Splitting Conic Solver (c) Brendan O'Donoghue, Stanford University, 2012 ------------------------------------------------------------------ problem: variables n: 5, constraints m: 6 cones: s: psd vars: 6, ssize: 1 settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07 alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1 max_iters: 100000, normalize: 1, rho_x: 1.00e-06 acceleration_lookback: 10, acceleration_interval: 10 compiled with openmp parallelization enabled lin-sys: sparse-direct-amd-qdldl nnz(A): 6, nnz(P): 0 ------------------------------------------------------------------ iter | pri res | dua res | gap | obj | scale | time (s) ------------------------------------------------------------------ 0| 1.97e+02 1.05e+01 4.39e+03 -2.20e+03 1.00e-01 1.51e-04 50| 6.29e-06 1.81e-06 2.22e-05 -1.11e-05 1.00e-01 2.99e-04 ------------------------------------------------------------------ status: solved timings: total: 3.00e-04s = setup: 4.59e-05s + solve: 2.54e-04s lin-sys: 1.37e-05s, cones: 1.30e-04s, accel: 3.29e-06s ------------------------------------------------------------------ objective = -0.000011 ------------------------------------------------------------------
This time, SCS reports 5
variables and 6
constraints.
The difference comes from the fact that, when designing the JuMP interface of
SCS, it was decided that the model would be read in the image form.
SCS therefore declares that it only supports free variables, represented in
JuMP as variables in MOI.Reals
and affine semidefinite constraints,
represented in JuMP as
MOI.VectorAffineFunction
-in-MOI.PositiveSemidefiniteConeTriangle
constraints.
On the other hand, SumOfSquares gave the model in kernel form so the
positive semidefinite (PSD) variables were reformulated as free variables
constrained to be PSD using an affine PSD constraints.
This transformation is done transparently without warning but it can be
inspected using print_active_bridges
.
As shown below, we can see
Unsupported variable: MOI.PositiveSemidefiniteConeTriangle
and
adding as constraint
indicating that PSD variables are not supported and they are added as free
variables.
Then we have Unsupported constraint: MOI.VectorOfVariables-in-MOI.PositiveSemidefiniteConeTriangle
indicating that SCS does not support constraining variables in the PSD cone
so it will just convert it into affine expressions in the PSD cone.
Of course, this is equivalent but it means that SCS will not exploit this
particular structure of the problem hence solving might be less efficient.
print_active_bridges(model_scs)
* Supported objective: MOI.ScalarAffineFunction{Float64} * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}} | bridged by: | SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, FullSpace, Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}, Vector{Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}}}}, Union{SumOfSquares.EmptyCone, SumOfSquares.PositiveSemidefinite2x2ConeTriangle, MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeTriangle, MonomialBasis, MonomialBasis, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}} | may introduce: | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-PolyJuMP.ZeroPolynomialSet{FullSpace, MonomialBasis, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}} | | bridged by: | | PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}} | | may introduce: | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Zeros | * Unsupported variable: SumOfSquares.EmptyCone | | adding as constraint: | | * Supported variable: MOI.Reals | | * Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.EmptyCone | | | bridged by: | | | SumOfSquares.Bridges.Constraint.EmptyBridge{Float64, MOI.VectorOfVariables} | | | may introduce: | * Unsupported variable: MOI.Nonnegatives | | adding as constraint: | | * Supported variable: MOI.Reals | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.Nonnegatives | | | bridged by: | | | MOIB.Constraint.FunctionConversionBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables, MOI.Nonnegatives} | | | may introduce: | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Nonnegatives | * Unsupported variable: SumOfSquares.PositiveSemidefinite2x2ConeTriangle | | bridged by: | | SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge{Float64} | | may introduce: | | * Unsupported variable: MOI.RotatedSecondOrderCone | | | adding as constraint: | | | * Supported variable: MOI.Reals | | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.RotatedSecondOrderCone | | | | bridged by: | | | | MOIB.Constraint.RSOCtoSOCBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables} | | | | may introduce: | | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.SecondOrderCone | * Unsupported variable: MOI.PositiveSemidefiniteConeTriangle | | adding as constraint: | | * Supported variable: MOI.Reals | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.PositiveSemidefiniteConeTriangle | | | bridged by: | | | MOIB.Constraint.SetDotScalingBridge{Float64, MOI.PositiveSemidefiniteConeTriangle, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables} | | | may introduce: | | | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle} | | | | bridged by: | | | | SCS.ScaledPSDConeBridge{Float64, MOI.VectorAffineFunction{Float64}} | | | | may introduce: | | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-SCS.ScaledPSDCone
With the dual version, we can see that variables in the PSD cone are supported directly hence we don't need that extra conversion.
print_active_bridges(model_dual_scs)
* Supported objective: MOI.ScalarAffineFunction{Float64} * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}} | bridged by: | SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, FullSpace, Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}, Vector{Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}}}}, Union{SumOfSquares.EmptyCone, SumOfSquares.PositiveSemidefinite2x2ConeTriangle, MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeTriangle, MonomialBasis, MonomialBasis, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}} | may introduce: | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-PolyJuMP.ZeroPolynomialSet{FullSpace, MonomialBasis, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}} | | bridged by: | | PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}} | | may introduce: | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Zeros | * Unsupported variable: SumOfSquares.EmptyCone | | adding as constraint: | | * Supported variable: MOI.Reals | | * Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.EmptyCone | | | bridged by: | | | SumOfSquares.Bridges.Constraint.EmptyBridge{Float64, MOI.VectorOfVariables} | | | may introduce: | * Supported variable: MOI.Nonnegatives | * Unsupported variable: SumOfSquares.PositiveSemidefinite2x2ConeTriangle | | bridged by: | | SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge{Float64} | | may introduce: | | * Supported variable: MOI.RotatedSecondOrderCone | * Supported variable: MOI.PositiveSemidefiniteConeTriangle
Consider a polynomial
$$
p(x) = \sum_{\alpha} p_\alpha x^\alpha,
$$
a vector of monomials b(x)
and the set
$$
\mathcal{A}_\alpha = \{\,(\beta, \gamma) \in b(x)^2 \mid x^\beta x^\gamma = x^\alpha\,\}
$$
The constraint encoding the existence of a PSD matrix Q
such that p(x) = b(x)' * Q * b(x)
can be written in standard conic form as follows:
$$
\begin{aligned}
\langle \sum_{(\beta, \gamma) \in \mathcal{A}_\alpha} e_\beta e_\gamma^\top, Q \rangle & = p_\alpha, \quad\forall \alpha\\
Q & \succeq 0
\end{aligned}
$$
Given an arbitrary choice of elements in each set $\mathcal{A}_\alpha$:
$(\beta_\alpha, \gamma_\alpha) \in \mathcal{A}_\alpha$.
It can also equivalently be written in the geometric conic form as follows:
$$
\begin{aligned}
p_\alpha e_{\beta_\alpha} e_{\gamma_\alpha}^\top +
\sum_{(\beta, \gamma) \in \mathcal{A}_\alpha \setminus (\beta_\alpha, \gamma_\alpha)}
y_{\beta,\gamma} (e_\beta e_\gamma - e_{\beta_\alpha} e_{\gamma_\alpha}^\top)^\top
& \succeq 0\\
y_{\beta,\gamma} & \text{ free}
\end{aligned}
$$
Let's study the evolution of the dimensions m
and n
of the semidefinite
program in two extreme examples and then try to extrapolate from these.
Suppose p
is a univariate polynomial of degree $2d$.
Then n
will be equal to d(d + 1)/2
for both the standard and geometric conic forms.
On the other hand, m
will be equal to 2d + 1
for the standard conic form and
d(d + 1) / 2 - (2d + 1)
for the geometric form case.
So m
grows linearly for the kernel form but quadratically for the image form!
Suppose p
is a quadratic form of d
variables.
Then n
will be equal to d
for both the standard and geometric conic forms.
On the other hand, m
will be equal to d(d + 1)/2
for the standard conic form and
0
for the geometric form case.
So m
grows quadratically for the kernel form but is zero for the image form!
In general, if $s_d$ is the dimension of the space of polynomials of degree d
then
$m = s_{2d}$ for the kernel form and $m = s_{d}(s_{d} + 1)/2 - s_{2d}$ for the image form.
As a rule of thumb, the kernel form will have a smaller m
if p
has a low number of variables
and low degree and vice versa.
Of course, you can always try with and without Dualization and see which one works best.
This notebook was generated using Literate.jl.